home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_13_11 / phillips / hurst.c < prev   
Text File  |  1995-01-18  |  8KB  |  275 lines
  1.  
  2.      /*******************************************
  3.      *
  4.      *   hurst(..
  5.      *
  6.      *   This routine performs the Hurst 
  7.      *   operation as described in "The Image 
  8.      *   Processing Handbook" by John C. Russ
  9.      *   CRC Press 1992.
  10.      *
  11.      *******************************************/
  12.  
  13. hurst(in_name, out_name, the_image, out_image,
  14.       il, ie, ll, le, size)
  15.    char   in_name[], out_name[];
  16.    int    il, ie, ll, le, size;
  17.    short  the_image[ROWS][COLS],
  18.           out_image[ROWS][COLS];
  19.  
  20. {
  21.    float  x[8], y[8], sig[8];
  22.    float  aa, bb, siga, sigb, chi2, q;
  23.    int    ndata, mwt;
  24.  
  25.    int    a, b, count, i, j, k,
  26.           new_hi, new_low, length,
  27.           number, sd2, sd2p1, ss, width;
  28.    short  *elements, max, prange;
  29.    struct tiff_header_struct image_header;
  30.  
  31.       /**********************************************
  32.       *
  33.       *   Initialize the ln's of the distances.
  34.       *   Do this one time to save computations.
  35.       *
  36.       **********************************************/
  37.  
  38.    x[1] = 0.0;        /* ln(1)        */
  39.    x[2] = 0.34657359; /* ln(sqrt(2))  */
  40.    x[3] = 0.69314718; /* ln(2)        */
  41.    x[4] = 0.80471896; /* ln(sqrt(5))  */
  42.    x[5] = 1.03972077; /* ln(sqrt(8))  */
  43.    x[6] = 1.09861229; /* ln(3)        */
  44.    x[7] = 1.15129255; /* ln(sqrt(10)) */
  45.  
  46.    sig[1] = 1.0;
  47.    sig[2] = 1.0;
  48.    sig[3] = 1.0;
  49.    sig[4] = 1.0;
  50.    sig[5] = 1.0;
  51.    sig[6] = 1.0;
  52.    sig[7] = 1.0;
  53.  
  54.    sd2 = size/2;
  55.  
  56.       /**********************************
  57.       *
  58.       *   Create out file and read
  59.       *   input file.
  60.       *
  61.       ***********************************/
  62.  
  63.    create_file_if_needed(in_name, out_name, out_image);
  64.  
  65.    read_tiff_image(in_name, the_image, il, ie, ll, le);
  66.  
  67.    max = 255;
  68.    if(image_header.bits_per_pixel == 4){
  69.       max = 16;
  70.    }
  71.  
  72.       /***************************
  73.       *
  74.       *   Loop over image array
  75.       *
  76.       ****************************/
  77.  
  78.    printf("\n");
  79.    for(i=sd2; i<ROWS-sd2; i++){
  80.       if( (i%2) == 0) printf("%d ", i);
  81.       for(j=sd2; j<COLS-sd2; j++){
  82.  
  83.          for(k=1; k<=7; k++) y[k] = 0.0;
  84.  
  85.             /*************************************
  86.             *
  87.             *   Go through each pixel class, set
  88.             *   the elements array, sort it, get
  89.             *   the range, and take the ln of the
  90.             *   range.
  91.             *
  92.             *************************************/
  93.  
  94.             /* b pixel class */
  95.          number      = 4;
  96.          elements    = (short *)
  97.                        malloc(number * sizeof(short));
  98.          elements[0] = the_image[i-1][j];
  99.          elements[1] = the_image[i+1][j];
  100.          elements[2] = the_image[i][j-1];
  101.          elements[3] = the_image[i][j+1];
  102.          sort_elements(elements, &number);
  103.          prange = elements[number-1] - elements[0];
  104.          if(prange < 0)  prange = prange*(-1);
  105.          if(prange == 0) prange = 1;
  106.          y[1] = log(prange);
  107.  
  108.             /* c pixel class */
  109.          elements[0] = the_image[i-1][j-1];
  110.          elements[1] = the_image[i+1][j+1];
  111.          elements[2] = the_image[i+1][j-1];
  112.          elements[3] = the_image[i-1][j+1];
  113.          sort_elements(elements, &number);
  114.          prange = elements[number-1] - elements[0];
  115.          if(prange < 0)  prange = prange*(-1);
  116.          if(prange == 0) prange = 1;
  117.          y[2] = log(prange);
  118.  
  119.             /* d pixel class */
  120.          elements[0] = the_image[i-2][j];
  121.          elements[1] = the_image[i+2][j];
  122.          elements[2] = the_image[i][j-2];
  123.          elements[3] = the_image[i][j+2];
  124.          sort_elements(elements, &number);
  125.          prange = elements[number-1] - elements[0];
  126.          if(prange < 0)  prange = prange*(-1);
  127.          if(prange == 0) prange = 1;
  128.          y[3] = log(prange);
  129.  
  130.             /* f pixel class */
  131.          if(size == 5  ||  size == 7){
  132.          elements[0] = the_image[i-2][j-2];
  133.          elements[1] = the_image[i+2][j+2];
  134.          elements[2] = the_image[i+2][j-2];
  135.          elements[3] = the_image[i-2][j+2];
  136.          sort_elements(elements, &number);
  137.          prange = elements[number-1] - elements[0];
  138.          if(prange < 0)  prange = prange*(-1);
  139.          if(prange == 0) prange = 1;
  140.          y[5] = log(prange);
  141.          }  /* ends if size == 5 */
  142.  
  143.             /* g pixel class */
  144.          if(size == 7){
  145.          elements[0] = the_image[i-3][j];
  146.          elements[1] = the_image[i+3][j];
  147.          elements[2] = the_image[i][j-3];
  148.          elements[3] = the_image[i][j+3];
  149.          sort_elements(elements, &number);
  150.          prange = elements[number-1] - elements[0];
  151.          if(prange < 0)  prange = prange*(-1);
  152.          if(prange == 0) prange = 1;
  153.          y[6] = log(prange);
  154.          }  /* ends if size == 7 */
  155.  
  156.          free(elements);
  157.  
  158.             /* e pixel class */
  159.          if(size == 5  ||  size == 7){
  160.          number      = 8;
  161.          elements    = (short *)
  162.                        malloc(number * sizeof(short));
  163.          elements[0] = the_image[i-1][j-2];
  164.          elements[1] = the_image[i-2][j-1];
  165.          elements[2] = the_image[i-2][j+1];
  166.          elements[3] = the_image[i-1][j+2];
  167.          elements[4] = the_image[i+1][j+2];
  168.          elements[5] = the_image[i+2][j+1];
  169.          elements[6] = the_image[i+2][j-1];
  170.          elements[7] = the_image[i+1][j-2];
  171.          sort_elements(elements, &number);
  172.          prange = elements[number-1] - elements[0];
  173.          if(prange < 0)  prange = prange*(-1);
  174.          if(prange == 0) prange = 1;
  175.          y[4] = log(prange);
  176.          }  /* ends if size == 5 */
  177.  
  178.             /* h pixel class */
  179.          if(size == 7){
  180.          elements[0] = the_image[i-1][j-3];
  181.          elements[1] = the_image[i-3][j-1];
  182.          elements[2] = the_image[i-3][j+1];
  183.          elements[3] = the_image[i-1][j+3];
  184.          elements[4] = the_image[i+1][j+3];
  185.          elements[5] = the_image[i+3][j+1];
  186.          elements[6] = the_image[i+3][j-1];
  187.          elements[7] = the_image[i+1][j-3];
  188.          sort_elements(elements, &number);
  189.          prange = elements[number-1] - elements[0];
  190.          if(prange < 0)  prange = prange*(-1);
  191.          if(prange == 0) prange = 1;
  192.          y[7] = log(prange);
  193.          }  /* ends if size == 7 */
  194.  
  195.          free(elements);
  196.  
  197.             /*************************************
  198.             *
  199.             *   Call the fit routine to fit the
  200.             *   data to a straight line. y=mx+b
  201.             *   The answer you want is the slope
  202.             *   of the line.  That is returned
  203.             *   in the parameter bb.
  204.             *
  205.             *************************************/
  206.          ndata = size;
  207.          mwt   = 1;
  208.          fit(x, y, ndata, sig, mwt, &aa, &bb,
  209.              &siga, &sigb, &chi2, &q);
  210.  
  211.          out_image[i][j] = (short)(bb*64.0);
  212.          if(out_image[i][j] > max)
  213.             out_image[i][j] = max;
  214.          if(out_image[i][j] < 0)
  215.             out_image[i][j] = 0;
  216.  
  217.       }  /* ends loop over j */
  218.    }  /* ends loop over i */
  219.  
  220.  
  221.    fix_edges(out_image, sd2);
  222.  
  223.    write_array_into_tiff_image(out_name, out_image,
  224.                                il, ie, ll, le);
  225.  
  226. }  /* ends hurst */
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.     /***********************************************
  234.     *
  235.     *    sort_elements(...
  236.     *
  237.     *    This function performs a simple bubble
  238.     *    sort on the elements from the median
  239.     *    filter.
  240.     *
  241.     ***********************************************/
  242.  
  243. sort_elements(elements, count)
  244.    int   *count;
  245.    short elements[];
  246. {
  247.    int i, j;
  248.    j = *count;
  249.    while(j-- > 1){
  250.       for(i=0; i<j; i++){
  251.          if(elements[i] > elements[i+1])
  252.             swap(&elements[i], &elements[i+1]);
  253.       }
  254.    }
  255. }  /* ends sort_elements */
  256.  
  257.  
  258.  
  259.     /***********************************************
  260.     *
  261.     *    swap(...
  262.     *
  263.     *    This function swaps two shorts.
  264.     *
  265.     ***********************************************/
  266.  
  267. swap(a, b)
  268.    short *a, *b;
  269. {
  270.    short temp;
  271.    temp  = *a;
  272.    *a    = *b;
  273.    *b    = temp;
  274. }  /* ends swap */
  275.